Landscape features influencing gene flow and connectivity of an endangered passerine

Abstract Dispersal of individuals and gene flow are crucial aspects to maintain genetic diversity and viability of populations, especially in the case of threatened species. Landscape composition and structure may facilitate or limit individual movement within and among populations. We used a landscape genetics approach to assess the connectivity patterns of the threatened Dupont's lark (Chersophilus duponti subsp. duponti), considering their genetic patterns and the landscape features associated with its gene flow in Spain. We analysed the genetic relatedness based on 11 species‐specific polymorphic microsatellites on 416 Dupont's lark individuals sampled across peninsular Spain between 2017 and 2019, covering most of the European distribution of the species. To assess the relationship between the landscape composition and the species gene flow, we estimated genetic distance at the individual level (Dps). Next, we built a set of environmental surfaces from two time periods (years 1990 and 2018), based on factors such as land use and topography, influencing individuals' movement. We then obtained resistance surfaces from an optimization process on landscape variables. Landscape genetics analyses were done for single and composite surface models for each year separately. Our findings from both time points show that scatter or mosaic‐structured vegetation composed by low agricultural and tree cover and high presence of sclerophyllous shrubs favoured Dupont's lark dispersal, while dense and continuous tree cover, as well as areas of intensive agriculture, were limiting factors. Our results suggest the importance of steppe habitat patches for the species' establishment and dispersal. In addition, our results provide key information to develop conservation measures, including conserving and restoring steppe habitats as scattered and/or mosaic‐structured vegetation that could warrant the connectivity and persistence of Dupont's lark populations.


| INTRODUC TI ON
From an ecological perspective, dispersal influences the distribution and abundance of species and the dynamics and persistence of populations.Besides, from an evolutionary perspective, dispersal affects processes like local adaptation and speciation since it determines the patterns of gene flow among populations (Dieckman et al., 1999;Peninston et al., 2023).Habitat quality, composition and structure are important factors affecting the dispersal, establishment and survival of individuals and hence, populations (Sunny et al., 2022).Individual dispersal, genetic diversity and gene flow, together with effective reproduction, are crucial for the viability of species (Balkenhol et al., 2015;Wittische et al., 2019) and their potential to adapt to environmental changes (Blair et al., 2013;Connor et al., 2022).The environmental features and habitat heterogeneity that define a landscape greatly influence dispersal patterns in most species, because they either facilitate or restrict the movement of individuals (Holderegger & Wagner, 2006).In turn, landscape connectivity and individual dispersal ability are factors that influence population differentiation, yielding the two most common patterns.
One, isolation by distance (IBD; Wright, 1943), entails increasing differentiation among populations as geographic distance increases; and the other, isolation by resistance (IBR; McRae, 2006), where population differentiation correlates with environmental resistance to individual movement while considering all possible pathways connecting them (McRae & Beier, 2007).Habitat fragmentation and loss modify the structure and composition of the landscape, creating habitat patches of different sizes, degrees of connectivity and spatial variability of resource availability (e.g.food, shelter).Some likely ecological and evolutionary consequences of habitat fragmentation and loss are declines in population size and genetic diversity and higher genetic structure (Alcaide et al., 2009;Balkenhol et al., 2015;García et al., 2021;Latorre-Cardenas et al., 2021), leading to local extinctions (Caizergues et al., 2003), which are expected to be more severe in species with restricted dispersal and strict habitat requirements.
Genetic methods are commonly used to indirectly estimate migration rates and effective dispersal, which have the advantage of capturing information about the movement of genes (Fukuda et al., 2022;Koopman et al., 2007).Specifically, the field of landscape genetics focuses on evaluating the influence of landscape features (composition, structure and matrix quality) on microevolutionary processes (i.e. gene flow, genetic drift and selection), at the individual or population level (Balkenhol et al., 2015;Manel et al., 2003;Storfer et al., 2007).Landscape genetics combines landscape ecology, population genetics and spatial statistics to determine the effects of geographic distance and landscape attributes on genetic diversity, structure and gene flow (Baden et al., 2019;Flores-Manzanero et al., 2019;Latorre-Cardenas et al., 2021;Shirk et al., 2017;Storfer et al., 2007).It allows the evaluation of functional connectivity, that is, the association between the genetic distance of genetically differentiated populations and the landscape characteristics.To achieve this, resistance surfaces are frequently used, and based on a variety of landscape features like vegetation, topography and climatic and anthropogenic variables, among others, to evaluate their relationship with genetic distances (IBR; McRae, 2006).The use of resistance models not only provides a means to assess the impact of the environment on population genetics but also provides insights into landscape management (e.g.design of corridors, detection of barriers), with applications to deal with current and future landscape changes and the maintenance and conservation of species (Balkenhol et al., 2015;Borokini et al., 2021;Shirk et al., 2017).
These habitats are currently experiencing strong degradation (Morales & Traba, 2013;Werger & van Staalduinen, 2012;Zurdo et al., 2021), mainly due to landscape transformation caused by human activity (e.g.agriculture intensification, wind farms and extensive grazing abandonment; Gómez-Catasús et al., 2021;Reverter et al., 2021).The result of these activities is habitat fragmentation and loss, with an increasing number of small steppe patches, which significantly limits the availability of Dupont's lark suitable habitat.Indeed, such habitat deterioration is one of the main factors of its Iberian range contraction (ca.36% in the last two decades) and population decline (ca.2.3% per year; Reverter et al., 2023).Recent estimations consider less than 2300 males in the whole Spanish population, with an effective population size of around 600-1300 pairs (Reverter et al., 2023).Although the impact of habitat degradation on genetic diversity is less evident in birds in comparison with other vertebrate species because of their higher dispersal abilities (Alcaide et al., 2009;Burgess & Garrick, 2020;Canales-Delgadillo et al., 2012;Coster et al., 2019;Ferrer et al., 2016;Mulvaney et al., 2021), bird species with restricted dispersal and strict habitat requirements are prone to be genetically affected by habitat fragmentation and loss, often exhibiting significant population structure and low genetic variation (Caizergues et al., 2003;Jensen et al., 2019;Khimoun et al., 2017;Klinga et al., 2019).
Previous studies reported variable levels of genetic structure across the range of Dupont's lark, depending on the geographic scale of analysis and the informativeness of the markers used.
Mitochondrial analyses at a range-wide scale found that the Spanish clade encompasses six closely related haplotypes, most of them shared across the Iberian Peninsula (García et al., 2008); thus, no meaningful patterns of mitochondrial differentiation were detected.
Microsatellite analyses reported moderate but statistically significant genetic differentiation for the Spanish populations (Méndez, Tella, & Godoy, 2011), rejecting the hypothesis of one single panmictic population.Some of the genetic clusters identified were unique for certain Spanish regions while others were shared among them with high gene flow.Such genetic structure partially coincides with five ecogeographic regions that differ in topography, climate, landscape distribution and amount of suitable habitat within and among clusters (Laiolo & Tella, 2006a, 2006b).The fact that the Dupont's lark populations show both certain structure in genetic clusters and shared mitochondrial haplotypes among regions indicates the existence of historical gene flow, while coalescent-based models and departures from mutation-drift equilibrium support a recent reduction in gene flow attributable, at least in part, to the transformation of steppes in the last decades (Méndez, Tella, & Godoy, 2011).This was recently corroborated by a study monitoring temporal trends in population genetic diversity and structure over the last 5-7 generations (Bustillo-de la Rosa et al., 2022).Authors found signs of genetic erosion and temporal changes in gene flow, which was overall high although variable among regions, and highly dependent on landscape-mediated connectivity among populations.In fact, the landscape configuration can also act as a potential barrier in this species, as shown by a study of song dissimilarities (Laiolo & Tella, 2006b), whereas suitable habitat patches work as steppingstones connecting inhabited areas, i.e. habitat patches identified as suitable for the species establishment and movement among occupied patches at distances up to 20 km (García-Antón et al., 2021).
Despite the importance of landscape on fragmentation and isolation among populations shown by the aforementioned studies, its relationship with genetic connectivity and gene flow has not been explicitly assessed.Furthermore, the demographic and genetic negative effects associated with population isolation (e.g.population size decline, genetic bottlenecks, diversity loss, genetic drift) can have crucial ecological and evolutionary consequences on, among others, metapopulation dynamics and local adaptation, increasing extinction risk especially in low dispersal specialist species (Peninston et al., 2023;Wong et al., 2023).Indeed, isolated populations resulting from land-use changes are of major concern in evolutionary biology, conservation and human and animal health (Vázquez-Domínguez et al., 2024).Therefore, and considering the current habitat modification of natural steppes and the reported decline of the species' populations, it is urgent to evaluate Dupont's lark functional connectivity and determine the landscape features associated with its dispersal and survival.
The aim of this study was to assess the relationship between the landscape composition and structure and the species' gene flow, identifying the main environmental factors influencing connectivity patterns of this threatened species.To this end, we considered two time periods to account for the potential time lag between landscape modifications and its genetic consequences over time.
Based on Dupont's lark habitat preferences and restricted dispersal (Gómez-Catasús et al., 2016;Suárez, 2010), we predicted that connectivity would be facilitated by the presence of steppe patches, characterized mainly by pasture and small shrubs.In contrast, that habitat where vegetation is conspicuously absent, with dense forests, or with agricultural fields would restrict movements.Regarding the anthropogenic landscape modification, potential differences at the two times evaluated in the landscape variables identified as influencing gene flow could be expected.Our results can provide key information to aid in slowing the current process of habitat degradation, further isolation and potential extinction of the species.
Namely, for addressing new insights into conservation measures, including the identification of habitat patches as suitable environmental corridors among populations that could conserve and improve the genetic diversity, connectivity and persistence of Dupont's lark populations.Considering the Dupont's lark as a bioindicator of steppe areas, identifying the key landscape variables favouring gene flow and survival can provide general information regarding ecological management of related steppe species.

| Study system and sampling
Dupont's lark individuals were sampled during the breeding season across peninsular Spain between 2017 and 2019 (Figure 1, see Figure A1 in Appendix A for the main landscape features within the study area), considering the locations where its presence was previously registered during the II National Census of the species (2004-2007;Suárez, 2010).Captures were carried out using individual spring-traps baited with mealworms (Tenebrio molitor), and fieldwork was performed exclusively by a team specialized in Dupont's lark.
We also used Dupont's lark recordings played through speakers to attract individuals; thus, sampling was biased towards adult males as described by Laiolo et al. (2007).Every individual captured was metal-ringed for identification and to avoid resampling.UTM coordinates of every individual captured were recorded using a GPS.
We collected blood samples from the jugular or brachial vein using 0.5 mL insulin syringes and capillary tubes; samples were stored in pure ethanol and transported to the laboratory for genetic analyses.
All individuals were released at the point of capture, keeping handling time as short as possible to avoid potential stress.The protocol for capture, handling, blood-sampling and release was approved by the Ethical Committee for Animal Experiments of the Universidad Autónoma de Madrid (CEI80-1468-A229).

| Microsatellites amplification and genetic analyses
In the present study, we used the microsatellite loci genotypes obtained by Bustillo-de la Rosa et al. (2022) where the DNA extraction protocol and PCR conditions are explained.Briefly, DNA was extracted from blood samples which were genotyped at 12 polymorphic microsatellite loci isolated in C. duponti (Méndez, Prats, et al., 2011).One locus showed no polymorphism and was removed from further analyses (Bustillo-de la Rosa et al., 2022).Amplifications were performed in one multiplex panel using four dyes (FAM, PET, NED, VIC) in a 10 μL total volume containing: 25-50 ng of template DNA, 1× QIAGEN Multiplex PCR Master Mix and 0.1 μM of each forward and reverse primers.PCR products were processed using capillary electrophoresis with an ABI 3730 by the Unidad de Genómica, Universidad Complutense de Madrid, Spain (https://www.ucm.es/ gyp/ genomica).Electropherograms were scored using Geneious v.10.2.3 (https://www.genei ous.com).Negative controls (i.e.ddH 2 O) were included in all runs and 7% of the samples were re-amplified and re-scored to confirm our genotypes and to estimate the genotyping error rate, which was less than 1% (see Table B1 in Appendix B for further information about the genotyping error rate per locus).The original data set was evaluated with Genepop on the web (Raymond & Rousset, 1995) and results showed that no microsatellite loci exhibited evidence of significant deviation from Hardy-Weinberg or linkage disequilibrium.
To avoid a potential bias from non-random sampling and overrepresentation of allele frequencies in our data due to the presence of closely related individuals (Waples & Anderson, 2017), we assessed relatedness between pairs of individuals with ML-Relate (Kalinowski et al., 2006).Less than 1% of pairs exhibited a first-order relationship (full-sibling, parent-offspring), thus all individuals were kept for analyses.Genetic population structure of Dupont's lark has been previously addressed with different genetic markers (Bustillo-de la Rosa et al., 2022;García et al., 2008;Méndez, Tella, & Godoy, 2011), showing that there are both unique and shared clusters in Spain with some degree of gene flow among them.Thus, landscape genetic analyses were performed at the individual level, calculating genetic distance between all pairs of individuals, increasing the likelihood to find genetic differences at finer scales and to detect novel barriers to gene flow (Mulvaney et al., 2021;Pavlacky Jr. et al., 2009).We estimated the proportion of shared alleles (Ps) with the function prop-Shared in adegenet and subtracted those values from one to obtain a genetic distance measure (Dps = 1 − Ps), which is an adequate metric for performing individual-based genetic distance estimates and useful for fine-scale landscape genetic analyses (Flores-Manzanero et al., 2019;Kimmig et al., 2020;McCluskey et al., 2022;Shirk et al., 2017).We tested for IBD with a Mantel test between genetic (Dps) and geographic distances; geographic (Euclidean) distance was estimated with GenAlEx v.6.5 (Peakall & Smouse, 2006, 2012).
Significance was obtained using 1000 permutations.
To incorporate explicit information about the species dispersal distances and address potential movement thresholds, we also performed Mantel correlograms based on mean dispersal distances (Storfer et al., 2007).We selected distance class boundaries of 5 and 20 km described for Dupont's lark resident movements and juvenile dispersal, respectively, which reflect sensible biological classes ( García-Antón et al., 2021).We employed an individual-based method using multi-locus spatial autocorrelation analyses in GenAlEx 6.5 (Banks & Peakall, 2012;Peakall & Smouse, 2006) to test the absence of spatial genetic structure.This method provides a measure of genetic correlation as a function of geographic distance.The method requires as inputs a geographic and a genetic distance matrix and provides an estimate of the autocorrelation coefficient (r) for each distance class.The null hypothesis (no spatial structure; r = 0) was obtained by 1000 random permutations of r values (random shuffling of all individuals among the geographic locations).If the calculated r values fall outside the 95% confidence interval bounding the null hypothesis, significant spatial genetic structure is inferred.We also defined the 95% confidence interval for each estimate of r by 999 bootstraps within each distance class (Peakall et al., 2003).

| Landscape data
To evaluate the influence of landscape features on gene flow, we considered environmental predictors that could influence Dupont's lark movement based on current knowledge of the species' biology and habitat preferences (Gómez-Catasús et al., 2016;Méndez et al., 2014;Suárez, 2010), and which were previously used to model the potential presence of the species (Table 1; García-Antón et al., 2019) and movements (García-Antón et al., 2021).Anderson et al. (2010) suggested that the grain size (i.e.pixel size) of landscape data should be based on the distance of resident movements described for the study species.An estimate of 5 km has been recognized as the distance threshold for resident movements in Dupont's lark and for subpopulation boundaries (García-Antón et al., 2021), as well as the limit of cultural similarity in song structure among individuals (Laiolo, 2008) and from capture-recapture data (Pérez-Granados & López-Iborra, 2015).Hence, we used a grain size of 5 km.
To build the environmental surfaces, we created a 5 × 5 km grid cell across Peninsular Spain and estimated the percentage cover of each landscape variable per 5 km pixel.Next, to avoid potential bias due to an overextended study area, we first created a buffer of 10 km around each sampling location (García-Antón et al., 2021).
Then, we defined the study area as the surface surrounding all buffer locations.Last, we clipped all the surfaces estimated across the Peninsular Spain using the study area as a mask (Figure 1) and keeping the same grid cells resolution.
Contemporary land use layers were obtained from CORINE data available for the year 2018 (https:// land.coper nicus.eu/ ).However, the effect of landscape features on genetic processes cannot always be detected due to the time lag between environmental changes and the species' genetic response (Epps & Keyghobadi, 2015;Pavlacky Jr. et al., 2009).Moreover, regarding anthropogenic landscape changes, the analysis of temporal effects is hampered by the fact that the alteration is often a continuous process over a considerable period, making it difficult to determine a precise date for such alterations.Since theory predicts that negative genetic effects would increase over time with increasing numbers of generations (Young & Clarke, 2000), we also evaluated the effect of landscape features on gene flow using the oldest available landscape data set (CORINE land cover from 1990).Accordingly, it will allow us to assess whether the effect of landscape variables on gene flow is persistent over time (i.e.landscape variables affecting gene flow in both time points).
For this, we reclassified the CORINE original 1990 and 2018 land cover categories into 12 different surface, selecting the variables we considered could limit or facilitate Dupont's lark dispersal and that were adequate to define both suitable breeding areas (as in García-Antón et al., 2019) and potential corridors (Table 1): urban, dryland farming (Agr_dry), other farming (Agr_other), farming+natural mosaic (Agr_nat), farming+forest mosaic (Agr_tree), natural grassland (Pasture), sclerophyllous shrubs (Shrub_scle), other shrubs (Shrub_other), forest, scarce vegetation, unvegetated and water bodies (Table 1).Sclerophyllous shrubs were differentiated from other shrubs as they have specifically been previously described as favourable for the species presence.Furthermore, 'Shrub_other' comprises moors and heathlands (Code 322 of CORINE) and transitional woodland/shrub (Code 324), which represent closed vegetation of great height, defined as unfavourable TA B L E 1 Land use and topographic variables (% cover and mean values in a 5 × 5 km cell respectively) considered for the landscape analyses.

Land use
Urban a Artificial surfaces (111,112,121 ,122,123,124,131,132,133, 141,142) Not The urban surface was used only for the 2018 models.
for the species.In contrast, sclerophyllous shrubs included bushy sclerophyllous vegetation mixed with some bare ground (Code 323 of CORINE), considered a priory as potentially favourable for Dupont's lark (Table 1).
We used vector layers to calculate cover (expressed as percentage) of every surface at each grid cell, based on which we built raster layers for every land use.For the topographical variables, we obtained data for elevation, its coefficient of variation and slope from the National Centre of Geographic Information with a 25 m resolution (DTM 25 m; https:// www.cnig.es/ ).Topographic values were incorporated to the 5 × 5 km grid with the bilinear interpolation method.All raster cells that fell outside the Spanish geographic border were kept as 'NoData' values and were not considered during the optimization process and analyses.Next, variables that were mostly absent across the study area, namely those that had values of 0 in 80% or more of the pixels, were removed (scarce vegetation, unvegetated, water bodies and the slope's coefficient of variation for both time points, whereas urban was removed only for the 1990 data).To avoid collinearity, we calculated the variance inflation among all variables using the VIF function in R.Those variables with high collinearity were individually removed until none showed VIF > 5 (Fox & Monette, 1992).Only 'other farming' for both 1990 and 2018 was eliminated, resulting in a final data set with 10 variables for 2018 and nine for 1990 (Figure 2; Table 1).

| Landscape genetic analyses
We performed landscape genetic analyses for both data sets (1990 and 2018) separately.We applied the optimization framework developed by Peterman et al. (2014) to determine the resistance values of our environmental surfaces with no a priori assumptions about the scale and direction of the resistance relationship.The ResistanceGA (Peterman, 2018), testing both IBD and IBR scenarios (Peterman, 2014(Peterman, , 2018)).ResistanceGA explores the parameter space and maximizes the relationship between pairwise landscape distances (least-cost or resistance) and the genetic distance (in this case Dps), using a genetic optimization algorithm (GA, Scrucca, 2013).It also uses generalized linear mixed-effect models fitted using maximum likelihood population effects parameterization (MLPE) to account for the non-independence error associated with pairwise distances (Clarke et al., 2002).
The optimization process comprises several steps.was conducted twice to ensure convergence (Peterman, 2018).

| Landscape genetics
Optimization and bootstrap model selection for the single and composite surfaces from ResistanceGA analyses showed that landscape variables that performed better than geographic distance alone (i.e.IBD) differed between 1990 and 2018.Results for 2018 data set showed that the two best supported models (increment in AIC less than two among them) after the bootstrapping were dry farming (Agr_dry) and farming and forest mosaic surface (Agr_tree) (39.98% and 29.74% of the times respectively; Table 2).These two variables were defined by Inverse-Ricker and Ricker functions respectively.The following wellsupported model was Composite 5 (7.17%;ΔAIC = 2.5) which included the two variables: Agr_dry and Agr_tree, with a higher contribution of the former (68%) than the latter (32%) (Table 2, see Tables D1 and D2 in Appendix D for parameter estimates from mixed effects models).
Three other models had moderate support: Slope (6.83%), Agr_nat (6.32%) and Forest (5.55%) (3.2 < ΔAIC < 4.9), while the remaining models, comprising also the geographic distance (ΔAIC = 7.4) showed low levels of support (<2%; ΔAIC 5.7-28.1).Resistance results indicated that moderate extent (10%-40%) of dry arable land (Agr_dry) had the lowest level of resistance for Dupont's lark, with increasing resistance below and above these values (Figure 4a).Conversely, the agricultureforest mosaic (Agr_tree) showed that moderate levels (10%-40%) had the highest resistance (Figure 4b).Moreover, these two variables were significant predictors of genetic distance in the generalized linear mixed-effects models (see Tables D1 and D2 in Appendix D).Flat areas showed low resistance, with resistance increasing exponentially as the slope increased (Figure 4c); also, resistance peaked around 10%-20% Agr_nat cover (Figure 4d).Regarding Forest, areas with less than 60% forest cover showed little or no resistance, while areas with >80% were classified with the highest resistance (Figure 4e).All other single and composite resistance surfaces showed low or null contribution to gene flow (Table 2).
The optimized pattern observed for Agr_dry was the same as in 2018, assigning low resistance only to areas around 20%-30% cover (Figure 4f).Additionally, values associated with Pasture showed lowest resistance close to 20% cover, and high resistance at all other cover percentages (Figure 4g).Forest had also the same pattern in 1990 as in 2018 with low resistance below 60% (Figure 4h).As observed for 2018, lowest resistance values were identified where Agr_tree was absent (Figure 4i), but with a steep increase in resistance as farming and forest mosaic surface rises.Finally, lower resistance values were found with both absence and >40% sclerophyllous shrub surface (Shrub_scle) (Figure 4j), whereas the other shrub surface (Shrub_other) exhibited increasing resistance with increasing shrub cover (Figure 4k).Considering the percentage top ranked models, the two time periods showed consistent resistance patterns for Agr_dry and Forest, but with contrasting effect regarding pasture surfaces.

| DISCUSS ION
By using a landscape genetics approach, we were able to discern connectivity patterns of Dupont's lark in Spain with fine detail and at the individual level.Our findings demonstrate how the landscape matrix and composition influence Dupont's lark genetic patterns, identifying how steppe patches (i.e.flat areas covered by low shrubs and pasture) are favourable for its effective dispersal.In contrast, the presence of forest areas, high vegetation and/or intensive agriculture seem to limit gene flow.
The spatial pattern observed where individuals tended to be genetically more similar within a distance of approximately 20-30 km than those further apart strongly supports the proposed steppingstone model of connectivity for Dupont's lark (García-Antón et al., 2021;Méndez et al., 2014) instead of more general models based on long dispersal distances (e.g.Coster et al., 2019).This pattern agrees with the genetic structure observed in Spain, where some genetic clusters are not entirely isolated (Méndez, Tella, & Godoy, 2011).Similar patterns have been identified for other bird species mainly because of their high dispersal capacities (Kleinhans & Willows-Munro, 2019;Koopman et al., 2007;Kunz et al., 2022).
To directly evaluate the landscape influence on the stepping-stone connectivity, future studies should consider a sampling design that specifically accounts for the habitat features of patches potentially serving as stepping stones (Duforet-Frebourg & Slatkin, 2016).
Landscape composition and structure changed through time, as can be seen in the overall land use differences (surfaces) between the two time points we evaluated (Figure 2).Such variation highlights the importance of considering different temporal data sets, when available, in landscape genetics to assess individual genetic responses in relation to the landscape changes during the time lag considered (Epps & Keyghobadi, 2015;Landguth et al., 2010;Manel & Holderegger, 2013).Some caution is needed when interpreting the influence of land use surfaces on gene flow, as our model results could be influenced by the differences in the CORINE classification methods.CORINE landcover categories of 1990 were generated from a human photo-interpretation of satellite images and orthophotos at 1:100,000 scale.In contrast, land use categories for 2018 were generated at a 1:25,000 scale with a generalization of SIOSE 2014, which is an Information System of Land Occupation in Spain, included in the Nation Plan for Territorial Observation.The latter provided higher accuracy in terms of land use polygons delimitation for the whole country (García-Álvarez & Olmedo, 2022;Martínez-Fernández et al., 2019).Nonetheless, we are confident of the pattern observed given the concordance with Dupont's lark distribution and life history.much of the species' range (see Figure E1 in Appendix E for resistance maps).Interestingly, the low resistance values at moderate agriculture cover could be related to fallow lands, which have been described previously as potentially beneficial for Dupont's lark dispersal during the non-breeding season (Suárez et al., 2006).
Moreover, García-Antón et al. (2019) found that habitat patches where dry farming cover was higher than 40% could be unsuitable habitat for this species.2.
High resistance costs for forested areas were also identified with both temporal data sets, indicating that gene flow is strongly limited on highly forested areas and by high-vegetation cover (Agr_tree and Shrub_other).This is consistent with the ecology and habitat preferences of this species (Gómez-Catasús et al., 2016;Suárez, 2010).
Interestingly, the 1990 surfaces were important to identify areas with low pasture cover, namely areas with ≤30% of natural grassland cover mixed with shrub vegetation, together with high cover of sclerophyllous shrubland, as landscape variables facilitating individual movement.This confirms our prediction that steppe patches would be a key habitat supporting gene flow.Therefore, considering the best models obtained for both time periods, landscape features better explaining Dupont's lark dispersal are consistent with those environmental variables described as key factors for its occurrence.Open landscapes with scattered and/or mosaic structured vegetation, including low agricultural surface and dispersed trees combined with high sclerophyllous shrub cover may ease movement.In contrast, areas of intensive agriculture and dense and continuous tree cover function as barriers.
This pattern was also showed by Pavlacky Jr. et al. (2009) for the Australian logrunner (Orthonyx temminckii), where habitat preferences determining its occurrence were also important for its dispersal.
We acknowledge that our results have potential limitations, first because we only analysed adult males, and this species might exhibit sex and/or age bias in dispersal, whereas association with specific landscape features could differ between sexes (Garza et al., 2005).
Secondly, other environmental variables not considered in this study could also influence Dupont's lark gene flow, like habitat quality (i.e. resource availability), patch size, behavioural traits (e.g.idiosyncratic low gene flow due to its sedentary behaviour) or demographic patterns (e.g.lower genetic differentiation and higher genetic diversity for larger population sizes; Kimmig et al., 2020;Méndez et al., 2014).
Nonetheless, we were able to identify patterns of landscape heterogeneity influencing Dupont's lark gene flow, and thus, its potential impact on the species dispersal and persistence.
In all, it is imperative to guide management efforts to conserve and restore steppe landscape features, characterized by scat-

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors have no relevant financial or non-financial interests to disclose.

DATA AVA I L A B I L I T Y S TAT E M E N T
I confirm that the dataset used for this study is available in a Dryad repository: https:// datad ryad.org/ stash/ share/ rjij3 9tgv_ wT-bNGUN 9xTHG so2T-nJHrM FmtYt hZlCA .

F
I G U R E 1 Map showing the sampling localities of Dupont's lark males (black dots) in Spain (N = 416) during the breeding season (2017-2019).The dashed line depicts a 10 km buffer around each sample and the square polygon comprises the entire study area.
resistance surfaces optimization framework(Peterman, 2018;Peterman et al., 2014) uses Ricker and Monomolecular equations to transform data, in combination with linear mixed-effects models and genetic algorithms, not making a priori assumptions regarding the scale and direction of the resistance relationship.The Monomolecular and Ricker are two exponential-based functions used for ecological modelling, differing in the curve shape of the relationship they are modelling.This curve is mainly determined by shape and magnitude parameters, which produce a saturating exponential (growth or decay) curve for the Monomolecular function, and a hump-shaped curve (skewed to right or left) for the Ricker function(Bolker, 2008).During the optimization process, the genetic algorithm searches all possible combinations of these parameters for transforming resistance surfaces, denoted by 'r'(Peterman, 2018;Peterman et al., 2014).We determined landscape resistance values that best explained pairwise genetic dissimilarity (genetic distance) between individuals with F I G U R E 2 Raster layers (5 × 5 km) showing the different land uses considered for the resistance models for the Dupont's lark in Spain, based on the CORINE data 1990 and 2018 (left and right respectively).The colour bar shows the percentage cover of each land use, the elevation (in metres) and the slope (in percentage) for the study area.
First, continuous landscape variables were individually transformed based on Ricker and Monomolecular functions (eight transformations in total per variable; Bolker, 2008) using the GA.prep function with the default parameters.Transformations are performed to maximize the relationship between pairwise landscape and genetic distances.Next, surfaces were individually optimized estimating pairwise effective distances among individuals using the commuteDistance function with gdistance in R (van Etten, 2014), which evaluates the length of a random roundtrip move that connects two nodes; it uses an eight-neighbour connection scheme to assess connectivity(Flores-Manzanero et al., 2019;van Etten, 2014).We used Akaike's information criterion value (AIC) during the optimization process as the objective function, which was determined from linear mixedeffects models (lmem) fitted by the MLPE parameterization using lme4 in R(Bates et al., 2014).We used Dps as the response variable, the pairwise resistance distances for each landscape surface as the dependent variable, and the identity of pairwise combination of individuals sampling localities as random effect.The best-optimized resistance surfaces were estimated based on the AIC corrected for small/finite populations (AICc) value(Burnham & Anderson, 2002).By default, ResistanceGA also calculates a distance-only (IBD) and a null model to test model performance.To control for the potential bias from the sampling design (type I error), and to test the robustness of our models, we conducted a pseudo bootstrap analysis (Resist.bootfunction in R).We randomly selected 75% of the samples without replacement and optimized surfaces were then fit to the subsample data.Using 10,000 iterations, we calculated for each model the average rank, average model weight, average marginal R 2 (proportion of the variance explained by the fixed factors) and the proportion (percentage) in which a surface was chosen as the best model.The resulting surfaces were then chosen to build composite surfaces and run a multivariable optimization.Optimization (single and composite surface models)

Finally, the best
replicates (overall best AICc) of single and composite models were selected to run a bootstrap model selection to obtain the average model rank, average model weight and the top ranked model (being the models with highest percentage the ones identified as best models).3 | RE SULTS 3.1 | Population genetics and fine-scale genetic structure A total of 416 Dupont's lark males were captured during this study (11% of the estimated number of Dupont's lark males in Spain; Traba et al., 2019; Figure 1), and their DNA successfully amplified for 11 microsatellite loci.The Mantel test showed a weak positive albeit significant correlation between genetic (Dps) and Euclidean geographic distances (r = .049;p = .001;see Figure C1 in Appendix C for the graph representation of the Mantel test).On the other hand, Mantel correlograms showed positive and significant correlations (r = .004-.019; p < .05) up to a distance of about 20-30 km between samples, where individuals tended to be genetically more similar than expected by chance (Figure 3).As the 5 and 20 km correlograms showed similar patterns, only the 20 km one is shown.F I G U R E 3 Correlogram showing the genetic correlation r (mean value and 95% confidence error bars by bootstrapping) as a function of 20 km distance bins (10 distance classes, 20 km each up to 200 km).Red lines are the 95% CI of a null hypothesis of random distribution of genotypes (U, upper; L, lower), calculated using 999 permutations.Asterisks indicate statistically significant values.Sample sizes (n; number of pairs compared) of each distance class are also indicated.

F
Single surface optimization response curves obtained with ResistanceGA for the relationship between landscape features and genetic distance (Dps) for Dupont's lark in Spain in 2018 (panel above) and 1990 (below).Each curve represents the resistance cost imposed by each landscape variable after the optimization procedure.Histograms represent the frequency of each resistance value.The transformation equation and the model support for each variable are described in Table , our findings support the threatened status of Dupont's lark and call for urgent conservation measures.In this study, we have shown for the first time how, despite the potential high dispersal of bird species (suggested by the low relatedness values obtained), Dupont's lark gene flow and future survival are tightly linked to the presence of shrub steppe areas (recognized as a unique and threatened habitat in Europe; Sainz Ollero, 2013), impacted by land use transformations.Previous studies showed that the viability of Dupont's lark overall metapopulation is maintained by the current large size of some populations and the connectivity among them (Bustillo-de la Rosa et al., 2022).However, the species is suffering a severe decline (Reverter et al., 2023) following a centripetal contraction due to the extinction of peripheral populations.The most peripheral populations are the most affected by habitat fluctuations, given they are smaller and more isolated (García-Antón et al., 2021; tered and mosaic-structured vegetation with predominance of low sclerophyllous shrublands, which provide both suitable breeding grounds and stepping-stones for this endangered bird and related steppe species.For instance, by limiting the expansion of intensive crops or afforestation, or recovering extensive livestock grazing, whose decline has a direct impact on steppe habitat because it favours forest recovery following secondary vegetation succession(Martínez-Valderrama et al., 2021;Reverter et al., 2019; Traba &   Pérez Granados, 2022).AUTH O R CO NTR I B UTI O N SDaniel Bustillo-de la Rosa: Data curation (lead); formal analysis (lead); investigation (equal); methodology (equal); writing -original draft (lead); writing -review and editing (equal).Adrián Barrero: Investigation (supporting); writing -review and editing (supporting).Juan Traba: Conceptualization (lead); funding acquisition (lead); writing -review and editing (supporting).Jesús T. García: Investigation (equal); writing -review and editing (supporting).Manuel B. Morales: Writing -review and editing (supporting).Ella Vázquez-Domínguez: Data curation (lead); formal analysis (lead); methodology (equal); supervision (lead); writingoriginal draft (lead); writing -review and editing (supporting).ACK N OWLED G EM ENTS We thank Julia Gómez-Catasús, Margarita Reverter, Julia Zurdo, Israel Hervás and all LIFE projects staff for their support in field and logistics work.We thank C. Pérez-Granados and G.M. López-Iborra for sharing of Dupont's lark samples.We also thank the Instituto de Investigación en Recursos Cinegéticos, IREC (CSIC-UCLM-JCCM), who provided the laboratory facilities for DNA analysis.This work was approved by the Local Ethical Committee for Animal Experiments of the Universidad Autónoma de Madrid (CEI80-1468-A229).All participants accepted to be part of the study.This study was supported by the European Commission (LIFE Ricotí project LIFE15-NAT-ES-000802 and LIFE Connect Ricotí project LIFE20-NAT-ES-000133) and the FPI-UAM fellowship from the Universidad Autónoma de Madrid that supported D. Bustillo-de la Rosa.This paper contributes to project REMEDINAL TE-CM (P2018/EMT4338).
Bootstrap model selection results for the single and composite surfaces that resulted from multisurface optimization on genetic distance (D PS ) for Dupont's lark in Spain.Note: Landscape variables evaluated from 2018 and 1990.The fitted transformation for individual surfaces optimization is also shown.Bestsupported models are indicated by higher 'top ranked %' (number of times the model chosen performed better than the rest considering bootstrapping).k=number of parameters fit in each model plus the intercept.Average weight = (the probability that a model is the best in the model set, averaged over 1000 bootstrap replicates).AICc = Akaike information criterion value generated from the MLPE mixed effects model corrected for small sample size; ΔAICc is the AICc value adjusted for the number of populations sampled and the number of parameters optimized.weight=weight of support for surface i, given the surfaces assessed.R 2 m and R 2 c are the marginal and conditional values, respectively, of the fitted MLPE model.Composite surfaces include variables as explained in TableD1of the Appendix D.